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Abstract. Chebyshev rational approximation can be a viable method to com- 
pute the exponential of matrices with eigenvalues in the vicinity of the negative 
real axis, and it was recently applied successfully to solving nuclear fuel burnup 
equations. Determining the partial fraction decomposition (PFD) coefficients 
of this approximation can be difficult and they have been provided (for ap- 
proximation orders 10 and 14) by Gallopoulos and Saad in "Efficient solution 
of parabolic equations by Krylov approximation methods", SIAM J. Sci. Stat. 
Comput., 13(1992). It was recently discovered that the order 14 coefficients 
contain errors and result in 10 2 times poorer accuracy than expected by theory. 
The purpose of this note is to provide the correct PFD coefficients for approx- 
imation orders 14 and 16 and to briefly discuss the approximation accuracy 
resulting from the erroneous coefficients. 



1. Chebyshev rational approximation 

This note concerns the computation of matrix exponential based on the Cheby- 
shev rational approximation method (abbreviated CRAM in [7]) on the negative real 
axis. In this approach, the matrix exponential e At is approximated by a rational 
matrix function f(At), where the rational function f{z) is chosen as the best ratio- 
nal approximation of the exponential function on the negative real axis K_. Let 
7Tfc fe denote the set of rational functions rk,k( x ) — Pk( x ) I 'Qk( x ) where p k and q k 
are polynomials of order k. The CRAM approximation of order k is defined as the 
unique rational function rk,k = Pk( x ) / Qk( x ) satisfying 

(1) sup \f k . k (x) - e x \= inf <^ sup \r k , k (x) - e x \\ . 

The asymptotic convergence of this approximation on the negative real axis is 
remarkably fast with the convergence rate 0(H~ k ), where H — 9.289 025 49... 
is called the Halphen constant [3]. It was recently discovered by Stahl and 
Schmelzer [11] that this convergence extends to compact subsets on the complex 
plane and also to Hankel contours in C \ K_ . The application of this approx- 
imation to computing the matrix exponential was originally made famous by 
Cody, Meinardus, and Varga in 1969 in the context of rational approximation 
of e~ x on [0, oo ) and it was recently resurfaced by Trefethen, Weideman, and 
Schmelzer [T^]. For self-adjoint and negative semi-definite matrices, the method is 
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guaranteed to yield an error bound in 2-norm that corresponds to the maximum 
error of the rational approximation on the negative real axis. This has also been 
the main context for scientific applications [31 [OJ [TO]. Recently, the method has 
also been successfully applied to non-self-adjoint matrices with eigenvalues near 
the negative real axis [3 [6] . These specific matrices arise from a reactor physics 
application, where the changes in nuclide concentrations due to radioactive decay 
and neutron-induced reactions are governed by a linear system x' = Ax known as 
the burnup equations. 

2. Partial fraction decomposition form 

The main difficulty in using CRAM for computing the matrix exponential is de- 
termining the coefficients of the rational function for a given k. In principle, the 
polynomial coefficients of pk and qk can be computed with Remez-type methods but 
this requires delicate algorithms combined with high-precision arithmetics. Fortu- 
nately, these coefficients have been computed to a high accuracy by Carpenter et al. 
for approximation orders k = 0, 1, . . . , 30 and they are provided in [JJ. In practical 
applications, however, it is usually advantageous to employ CRAM in the partial 
fraction decomposition (pfd) form. For simple poles, this composition takes the 
form 



(2) fk,k( z ) = a + 



k 
\ " 



where is the limit of the function f k ,k at infinity, and aj are the residues at the 
poles 9j\ 



(3) 



Since the coefficients of f^k are real, its poles form conjugate pairs, so the compu- 
tational cost can be reduced to half for a real variable x 



(4) f*,*(a;) = oo+2Re J] 



fc/2 

a 4 



x — i 

3=1 



and the matrix exponential solution may be approximated as 

/fc/2 

(5) x = apxp + 2 Re [ otj(Ab — 9jl)~ x xp 

for a real matrix A GM. nxn . 

Although the PFD coefficients can in principle be computed from the polynomial 
coefficients, the computation of the polynomial roots may be ill-conditioned and 
requires great care. The PFD coefficients for approximation orders 10 and 14 have 
been provided in [3] , and the given coefficients for k = 14 have been used in several 
applications including the matrix exponential computing package expokit [TO] and 
the reactor physics code Serpent 5 . However, in the latter context, it was recently 
observed that these reported coefficients contain errors and do not correspond to 
the true best approximation [6 ]. To illustrate this, Figure [JJ shows the error of 
order 14 approximation on the negative real axis computed using two different sets 
of coefficients: the partial fraction coefficients from 3 a , with the corresponding 



CORRECTION TO CHEBYSHEV APPROXIMATION COEFFICIENTS 



3 





(b) 

FIGURE 1. (a) Plot of e x — 7*14,14(3;) on the negative real axis with 
ru,i4 computed based on the partial fraction coefficients from [3], (b) 
Plot of e x — fi4,i4(x) based on the polynomial coefficients from pQ. The 
plots were computed using high-precision arithmetics with 32 digits. 



approximation denoted by fi4,i4, and the polynomial coefficients from [T], with the 
corresponding approximation denoted by 7^14^4. According to theory, a necessary 
and sufficient condition for the best approximation is that the corresponding error 
function equioscillates, i.e. there exists a set of points where it attains its maximum 
absolute value with alternating signs. Notice that the approximation computed with 
the coefficients from [3] does not exhibit this behavior and in addition results in a 
10 2 times poorer accuracy than expected by theory. 

After discovering the erroneous behavior induced by the coefficients from [3], 
partial fraction coefficients for approximation orders k — 14 and k = 16 were com- 
puted from the polynomial coefficients provided in [1] and subsequently reported 
in [6] . The computed pfd coefficients are repeated here in Tables Q] and [2] The 
computations were performed with Matlab's Symbolic Toolbox using high preci- 
sion arithmetics with 200 digits to ensure a sufficient accuracy. In Tables [1] and [2] 
the coefficients have been rounded off to 20 digits. The coefficients in p] have been 
also given with 20 digits' accuracy, and based on our experience, the approximation 
order k = 16 is the highest for which this accuracy is sufficient for computing the 
pfd coefficients. For lower approximation orders, 1 < k < 13, the PFD coefficients 
can be accurately computed with the approximative Caratheodory-Fejer method 
and a Matlab script is provided for this purpose in [8]. 



4 



M. PUSA 



TABLE 1. Recomputed values for the partial fraction decomposition 
coefficients for CRAM approximation of order 14. 



Coefficient 


Real part 


Imaginary part 


8i 


-8.897 773 186 468 888 819 9 x 10° 


+1.663 098 261 990 208 530 4 x 10 1 


02 


-3.703 275 049 423 448 060 3 x 10° 


+1.365 637187148 326 8171 x 10 1 


83 


-0.208 758 638 250 130 125 1 x 10° 


+1.099 126 056 190 126 091 3 x 10 1 


04 


+3.993 369 710 578 568 519 4 x 10° 


+6.004 831642 235 037 317 8 x 10° 


05 


+5.089 345 060 580 624 506 6 x 10° 


+3.588 824 029 027006 510 2 x 10° 


0ti 


+5.623 142 572 745 977 124 8 x 10° 


+1.194069 046 343 966 9766 x 10° 


e 7 


+2.269 783 829 231 112 709 7 x 10° 


+8.461737 973 040 2214019 x 10° 


Oil 


-7.154 288 063 589 067 285 3 x 10~ 5 


+1.436104 334 954130 0111 x 10~ 4 


012 


+9.439 025 310 736 168 877 9 x 10~ 3 


-1.718 479 195 848 301 751 1 x 10~ 2 


a 3 


-3.763 600 387 822 696 871 7 x 10" 1 


+3.351834 702 945 010 4214 x lO" 1 


Ct4 


-2.349 823 209108 2701191 x 10 1 


-5.808 359129 714 207 400 4 x 10° 


0-5 


+4.693 327 448 883 129 304 7 x 10 1 


+4.564 364 976 882 776 0791 x 10 1 


ae 


-2.787 516194 014 564 646 8 x 10 1 


-1.021473 399 905 645 143 4 x 10 2 


a 7 


+4.807 112 098 832 508 890 7 x 10° 


-1.320 979 383 742 872 3881 x 10° 


ao 


+1.832 174 378 254 041275 1 x 10" 14 


+0.000 000 000 000 000 000 x 10° 



TABLE 2. Computed values for the partial fraction decomposition co- 
efficients for CRAM approximation of order 16. 



Coefficient 


Real part 


Imaginary part 


01 


-1.084 391707869 698 802 6 x 10 1 


+1.927 744 616 718165 228 4 x 10 1 


02 


-5.264 971343 442 646 889 5 x 10° 


+1.622 022147 316 792 730 5 x 10 1 


03 


+5.948152 268 951 177480 8 x 10° 


+3.587 457 362 018 322 282 9 x 10° 


04 


+3.509103 608 414 918 0974 x 10° 


+8.436198 985 884 375 082 6 x 10° 


05 


+6.416 177 699 099 434 192 3 x 10° 


+1.194122 393 370138 6874 x 10° 


0fi 


+1.419 375 897 185 665 978 6 x 10° 


+1.092 536 348 449 672 258 5 x 10 1 


07 


+4.993 174 737 717 996 3991 x 10° 


+5.996 881713 603 942 226 x 10° 


0s 


-1.413 928 462 488 886 2114 x 10° 


+1.349 772 569 889 274 538 9 x 10 1 


a\ 


-5.090 152 186 522 491 565 x 10~ 7 


-2.422 001765 285 228 797 x 10~ 5 


a2 


+2.115 174 218 246 603 090 7 x 10~ 4 


+4.389 296 964 738 067 3918 x 10~ 3 


«3 


+1.133 977517848 393 0527 x 10 2 


+1.019 472170 421585 645 x 10 2 


a 4 


+1.505 958 527 002 346 752 8 x 10 1 


-5.751405 277 642181997 9 x 10° 


0-5 


-6.450 087802 553 964 659 5 x 10 1 


-2.245 944 076 265 209 605 6 x 10 2 


«6 


-1.479 300 711 355 799 971 8 x 10° 


+1.768 658 832 378 293 790 6 x 10° 


a 7 


-6.251 839 246 320 791 889 2 x 10 1 


-1.119 039109 428 322 848 x 10 1 


as 


+4.102313683 5410021273 x 10~ 2 


-1.574 346 617 345 546 8191 x lO" 1 


ao 


+2.124 853 710 495 223 748 8 x lO" 16 


+0.000 000 000 000 000 000 x 10° 
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logio \?u,u( z ) - r u.u{ z ) 




Real axis 



FIGURE 2. Plot of log 10 |fi4,i4(z) — ^4,14(2)1 in the complex plane. 
?14,14 was computed based on the partial fraction coefficients from Ta- 
ble [1] and ri4 ; i4 was formed by truncating these coefficients to 6 sig- 
nificant digits. The poles of ri.4,14 have been marked to the plot with 
asterisks. 

3. Analysis of inaccurate pfd coefficients for k = 14 

To analyze the effect of inaccurate PFD coefficients denoted by {#_,} and {5j}, 
let r denote the corresponding rational approximation. The error caused by the 
inaccuracies in the PFD coefficients may be estimated 

k 11 ^ 

(6) \f k ,k( z ) ~ rk,k( z )\ S \a ~a \+Y^ i _\ 12 \°o ~ 9 j\ + i _ fl \ a i ~ S jl > 

indicating that the error is the greatest in the vicinity of the poles. It can also be 
seen from Eq. ^ that the inaccuracy related to the poles has a greater impact near 
the poles, whereas the error related to the residues should begin to dominate the 
total error farther away from the poles. By comparing the old and the recomputed 
PFD coefficients for k — 14, it can be seen that the poles all agree to about 6 digits 
whereas the residues agree to about 5 digits. El The most dramatic discrepancy 
occurs for the coefficient ao for which the significands agree to 5 digits but the 
exponent value given in [3] is —12, although the correct value is —14. 

On the grounds of Eq. ([6]), it can be estimated that coefficients with 6 correct 
digits should produce a rational function whose deviation from f 1 4 1 4(z) is at most 
of the order of 10 -3 on the negative real axis. Figure [2] shows the difference between 
^14 14 (z) an d a rational function r 14,14(2) that was formed by truncating the PFD 
coefficients of f 14,14 to 6 significant digits. Interestingly, as can be seen from Fig.QJi, 
the approximation #^4 14 (x), corresponding to the PFD coefficients from [3], yields a 
significantly better accuracy of order 10 -12 than is expected based on the accuracy 
of the coefficients alone. 



^Notice that the PFD coefficients in [3] are given for the rational approximation of e 1 on [0, 00) 
and that they have been multiplied by a factor of two making Eq. (37) in [3] equivalent to Eq. |(5jl. 
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To investigate the matter further, let us now take the poles {Oj} reported in [3] 
as a starting point for constructing a rational approximation of order 14. The poles 
{Oj} define a polynomial 

14 

(7) qii{x) = \{(x-e ] ) 

whose values agree to about 6 digits with the values of the correct polynomial qu (x) 
on the negative real axis. The residues at the poles {Oj} cannot be computed in a 
fully consistent manner, since the poles do not correspond to the true zeros of $14. 
However, two alternative approaches for computing the residues can be considered. 
One possibility is to use the correct rational function f 14^4 and Eq. ([3]) to compute 
the residues, but this is inconsistent as Eq. ([3]) only holds at the true poles. Another 
option is to define a new rational function using qu as the denominator and the 
correct polynomial p±4 as the numerator, after which the residues can be computed 
exactly using symbolic arithmetics. With both of these approaches we obtain a 
rational approximation, whose accuracy is of the order of 10 -6 on the negative 
real axis. It is also worth mentioning that forming the rational function based on 
the poles {Oj} and the correct residues {otj} from Table Q] yields an approximation 
whose accuracy is of the order of 10 -7 on the negative real axis. 

The article [3J by Gallopoulos and Saad does not indicate, how the reported PFD 
coefficients were computed, but based on the observations regarding the accuracy 
of the resulting approximation, it is evident that the values given for the residues 
somehow compensate for the inaccuracies in the poles {Oj} and it seems likely that 
they have been optimized to minimize the deviation from fi^n on the negative 
real axis. In fact, using the poles {Oj} and standard least squares optimization in 
Matlab with 10 7 points chosen from the interval [— 10 3 , — 10~ 10 ], we were able to 
produce residues yielding only a slightly worse accuracy of order 10 -11 . In any case, 
it should be noted that optimizing the residues properly in the Chebyshev sense 
would essentially form a problem of comparable difficulty as the original problem 
of determining 714,14. 
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